Pharmacokinetics

• Compartmental models
• Models with analytical solutions

Dosing regimens

Dosing regimens

Population simulation

PKPD

. # A tibble: 1 x 1
.   Median_CFB
.        <dbl>
. 1      -14.4

Sensitivity analysis

Mechanistic models

About mrgsolve

• R package for simulation from ODE-based models
  • Free, OpenSource, GitHub, CRAN
• Language
  • Models written in C++ inside model specification format
  • General purpose solver: ODEPACK / DLSODA (FORTRAN)
    • Automatically detect and switch between non-stiff (Adams) and stiff (BDF) methods for solving the differential equations
  • Simulation workflow in R
• Hierarchical (population) simulation
  • ID, \(\eta\), \(\varepsilon\)
• Integrated PK functionaility
  • Bolus, infusion, F, ALAG, SS etc, handled under the hood
  • 1- and 2-cmt PK models in closed-form
• Extensible using R, C++, Rcpp, boost, RcppArmadillo
• R is it’s natural habitat

mrgsolve started as QSP modeling tool

• Motivation: large bone/mineral homeostatsis model (CaBone)
• History using
  • Berkeley Madonna
  • WinBUGS
  • NONMEM (attempted)
• 2010: write R front end to deSolve
• 2012: write C++ interface to DLSODA
• Develop dosing / event capability
• More recently, expose functionality provided by
  • Rcpp - vectors, matrices, functions, environments, random numbers
  • boost - numerical tools in C++
  • users’ own C++ code (functions, data structures, classes)
• Translator from SBML to mrgsolve using R bindings to libSBML

Orientation

• https://CRAN.R-project.org/package=mrgsolve

• GitHub site: https://github.com/metrumresearchgroup/mrgsolve

• Issues and questions: https://github.com/metrumresearchgroup/mrgsolve/issues

• mrgsolve website: https://mrgsolve.github.io

• User Guide: https://mrgsolve.github.io/user_guide

• Blog: https://mrgsolve.github.io/blog

• Vignettes: https://mrgsolve.github.io/vignettes

What we will cover today

1. Three basic workflows
2. Loading the model into R
3. Event objects
4. Data sets

Emphasis is on getting you running your own simulations today.

A basic simulation with mrgsolve

mod %>% ev(amt = 100, ii = 24, addl = 3) %>% mrgsim() %>% plot()

The first / simplest workflow.

A basic simulation with mrgsolve

mod %>% ev(amt = 100, ii = 24, addl = 3) %>% mrgsim() %>% plot()

• mod: the model object
  • Ok … where did that come from?
• ev(amt = 100, …) : the intervention
  • An event in this example
• mrgsim(): actually do the simulation
• plot(): do something with the simulation
  • plot, mutate, as_data_frame etc …

• model %>% intervention %>% simulate %>% take-a-look

Why do we use %>% ?

What happens first in this operation?

mean(sqrt(seq(4)))

Pipelines

mean(sqrt(seq(4)))

4  %>% seq() %>% sqrt() %>% mean()

Better.

4  %>% 
  seq(.) %>% 
  sqrt(.) %>% 
  mean(., na.rm = TRUE)

mod %>% some_intervention() %>% simulate() %>% post_process()

The model object

model %>% ...

• I (almost) always call the model object mod in the documention / examples

• All the information about the model we need to know to run the simulation

• Distinct from the intervention, the population, the summary

• But the model does know about output time, random effects

Take a look: overview

mod

----------------  source: pk1.cpp  ----------------

  project: /Users/kyleb/Rli...gsolve/models
  shared object: pk1-so-d2b54722a2a3 

  time:          start: 0 end: 192 delta: 0.2
                 add: <none>

  compartments:  EV CENT [2]
  parameters:    CL V KA [3]
  captures:      CP [1]
  omega:         0x0 
  sigma:         0x0 

  solver:        atol: 1e-08 rtol: 1e-08 maxsteps: 5000

Take a look: parameters (really important)

param(mod)

 Model parameters (N=3):
 name value . name value
 CL   1     | V    20   
 KA   1     | .    .    

• Parameters get a name
• Names and number of parameters gets fixed at compile time
• All parameters have a value
• Value can be modified after compile time

Any advanced simulation will need to manipulate the model parameters in some way or another

• Can be CL, V2, Q etc
• Fixed effect parameter estimates (THETAs)
• Covariates (WT, EGFR etc)
• Flags
  • Is the dosing weight-based or flat? You could make a flag for that.

Take a look: compartments

init(mod)

 Model initial conditions (N=2):
 name       value . name     value
 CENT (2)   0     | EV (1)   0    

• Every compartment gets a name
• Every compartment gets an initial condition

Where did mod come from?

mod <- mread("simple", "../../model")

mod <- mread("<model-name>", "<project-directory>")

• By default mrgsolve looks for the code in the file
  • model-name.cpp in
  • project-directory
• mread demo

Internal model library

Quiz:

mod <- mread("<first-argument>", "<second-argument>")

Internal model library

Quiz:

mod <- mread("<first-argument>", "<second-argument>")

mod <- mread("<first-argument>", modlib())

modlib()

. [1] "/Users/kyleb/Rlibs/mrgsolve/models"

Internal model library

mod <- mread("effect", modlib())

mod

. 
. 
. ---------------  source: effect.cpp  ---------------
. 
.   project: /Users/kyleb/Rli...gsolve/models
.   shared object: effect-so-edd44d648855 
. 
.   time:          start: 0 end: 36 delta: 0.1
.                  add: <none>
. 
.   compartments:  GUT CENT PERIPH CE [4]
.   parameters:    VC KA K10 K12 K21 E0 EMAX EC50 KEO [9]
.   captures:      EFFECT CP [2]
.   omega:         0x0 
.   sigma:         0x0 
. 
.   solver:        atol: 1e-08 rtol: 1e-08 maxsteps: 5000

Try this now

mod <- mread("pbpk", modlib())

mod <- modlib("pbpk")

Oops

mod <- try(mread("simple"))

. Error : the model file simple.cpp does not exist.

mod

. [1] "Error : the model file simple.cpp does not exist.\n"
. attr(,"class")
. [1] "try-error"
. attr(,"condition")
. <simpleError: the model file simple.cpp does not exist.>

Your turn

Summary:

• Use mread to read, compile and load a model
• We can get information out of the model object that was returned
• This is the first step to get simulating (need a model first!)
• We will cover writing the model later on

• File name: exercises/warm_up.R

Roadmap

We use mread() to get a model loaded into R

Next, we’ll need

• An intervention to apply to the model
• A way to simulate from the model
• A way to deal with the output (like plot)

Event objects

• model %>% intervention %>% simulate %>% take-a-look

• Event object = quick / easy way to implement dose or other intervention

e <- ev(amt = 100) 

e

. Events:
.   time cmt amt evid
. 1    0   1 100    1

• Defaults: time, evid, cmt

What to include in ev(...)

• time event time
• cmt Event compartment
• amt Dose amount
• ii Inter-dose interval
• addl Additional doses to administer
• total Total number of doses to administer
• rate Infusion rate
• tinf Infusion duration
• ss Set to 1 to advance to steady state
• evid Event id
• ID Subject ID (use multiple ids - ID=1:10)

Interventions and corresponding evid

• Bolus dosing (evid 1, with rate==0)
• Zero order infusion (evid 1, with rate > 0)
• Other type event (evid 2)
  • This also forces solver reset
• Compartment reset (evid 3)
• Reset and dose (evid 4)
• Replace the amount in a specific compartment (evid 8)

demo

Three ways to invoke

Inline

mod %>% ev(amt = 100) %>% mrgsim()

Object via pipeline

e <- ev(amt = 100)

mod %>% ev(e) %>% mrgsim()

As argument

mod %>% mrgsim(events = e)

Event objects are just data frames

as.data.frame(e)

.   time cmt amt evid
. 1    0   1 100    1

• We will use a data_set later on for populations

• Event objects are convenient
  • Constructor
  • Operations

Simulate

model %>% intervention %>% simulate

e <- ev(amt = 100)

mod %>% ev(e) %>% mrgsim()

. Model:  effect 
. Dim:    362 x 8 
. Time:   0 to 36 
. ID:     1 
.     ID time    GUT   CENT PERIPH     CE EFFECT     CP
. 1:   1  0.0   0.00  0.000 0.0000 0.0000  157.0  0.000
. 2:   1  0.0 100.00  0.000 0.0000 0.0000  157.0  0.000
. 3:   1  0.1  91.21  8.443 0.1552 0.2225  155.7  3.460
. 4:   1  0.2  83.19 15.502 0.5817 0.8048  152.8  6.353
. 5:   1  0.3  75.88 21.354 1.2272 1.6382  149.6  8.752
. 6:   1  0.4  69.21 26.156 2.0463 2.6356  146.6 10.719
. 7:   1  0.5  63.13 30.046 3.0001 3.7278  144.1 12.314
. 8:   1  0.6  57.58 33.146 4.0553 4.8607  142.2 13.584

mrgsim(mod, events = e)

Simulate and plot

We need to know this to work the example; more details in a bit

mod %>% ev(e) %>% mrgsim() %>% plot()

Your turn

• Background

  • https://www.ncbi.nlm.nih.gov/pubmed/16988206
  • Meropenem is a carbapanem antibiotic
  • Like penicillins, efficacy related to time above MIC
  • The authors wanted to look at extended infusion durations

• File name: exercises/meropenem.R

Customizing output

What would you like to “fix” in this plot?

mod %>% ev(amt = 100) %>% mrgsim() %>% plot(CP~time)

Simulation time grid

• To form a sequence of times
  • start usually 0
  • end time of last observation
  • delta output time step
• Additional times to simulate
  • add ad-hoc numeric vector

stime(mod)

. [1]  0  6 12

Simulation end time

mod %>% ev(amt = 100) %>% mrgsim(end = 48) %>% plot(CP~time)

Simulation time step

mod %>% ev(amt = 100) %>% mrgsim(end = 48, delta = 0.1) %>% plot(CP~time)

Working with simulated output

model %>% intervention %>% simulate %>% take-a-look

Plot

mod <- mread_cache("effect", modlib())

mod %>% ev(amt = 100) %>% mrgsim() %>% plot()

Plot

mod %>% ev(amt = 100) %>% mrgsim() %>% 
  plot(CP + EFFECT~., col = "firebrick")

Pipeline to dplyr functions

mod %>% ev(amt = 100) %>% mrgsim() %>% mutate(arm = 1)

. # A tibble: 362 x 9
.       ID  time   GUT  CENT PERIPH    CE EFFECT    CP   arm
.    <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl> <dbl>
.  1     1   0     0    0     0     0       157   0        1
.  2     1   0   100    0     0     0       157   0        1
.  3     1   0.1  91.2  8.44  0.155 0.222   156.  3.46     1
.  4     1   0.2  83.2 15.5   0.582 0.805   153.  6.35     1
.  5     1   0.3  75.9 21.4   1.23  1.64    150.  8.75     1
.  6     1   0.4  69.2 26.2   2.05  2.64    147. 10.7      1
.  7     1   0.5  63.1 30.0   3.00  3.73    144. 12.3      1
.  8     1   0.6  57.6 33.1   4.06  4.86    142. 13.6      1
.  9     1   0.7  52.5 35.6   5.18  5.99    141. 14.6      1
. 10     1   0.8  47.9 37.4   6.36  7.09    139. 15.3      1
. # … with 352 more rows

If we didn’t modify the output

mod %>% mrgsim()

. Model:  effect 
. Dim:    361 x 8 
. Time:   0 to 36 
. ID:     1 
.     ID time GUT CENT PERIPH CE EFFECT CP
. 1:   1  0.0   0    0      0  0    157  0
. 2:   1  0.1   0    0      0  0    157  0
. 3:   1  0.2   0    0      0  0    157  0
. 4:   1  0.3   0    0      0  0    157  0
. 5:   1  0.4   0    0      0  0    157  0
. 6:   1  0.5   0    0      0  0    157  0
. 7:   1  0.6   0    0      0  0    157  0
. 8:   1  0.7   0    0      0  0    157  0

Other connections to dplyr

• mod %>% mrgsim() %>% mutate()
• mod %>% mrgsim() %>% filter()
• mod %>% mrgsim() %>% group_by()
• mod %>% mrgsim() %>% as_tibble()

Your turn

• Go back to the meropenem example
• Simulate 7 days of dosing 0.5 g Q 8 h over 0.5 hour
• Calcualte the fractional time above MIC on day 2 (MIC 2 or 4 mg/L)

• Demo: try 0.5 g Q 12 h over 0.5 hour for 1 day in normal and impaired renal function

• File name: exercises/meropenem.R

Updating the model object

Update

On the fly

mod %>% update(end = 72) %>% mrgsim()

Persistent

mod <- update(mod, end = 72)

mrgsim will call update for you (on the fly)

mod %>% mrgsim(end = 72)

What else can I update?

• Time
  • start, end, delta, add
• Parameters and compartment initial values
• Solver settings
  • atol, rtol
  • hmax, maxsteps, mxhnil, ixpr
• $OMEGA, $SIGMA
• tscale (rescale the output time)
• digits

mod %>% update(rtol = 1E-12) %>% ...

mod %>% mrgsim(rtol = 1E-12) %>% ...

Updating parameters

• Demo

mod <- param(mod, CL = 1.5)

Add a population element to the simulation

• model %>% intervention %>% simulate %>% take-a-look

• model %>% intervention %>% population %>% simulate %>% take-a-look

This is our second simulation workflow.

One population simulation with mrgsolve

mod %>%
ev(amt = 100, ii = 24, addl = 3) %>%
idata_set(pop) %>%
mrgsim() %>% plot()

idata_set takes in individual-level data

• ID - one per row
• Typically parameters in columns

head(pop, n = 3)

.   ID   CL   VC     KA KOUT IC50 FOO
. 1  1 1.05 47.8 0.8390 2.45 1.28   4
. 2  2 0.73 30.1 0.0684 2.51 1.84   6
. 3  3 2.82 23.8 0.1180 3.88 2.48   5

length(unique(pop$ID))

. [1] 10

• This tells mrgsolve to simulate 10 units or individuals

pop and mod are connected via parameters

head(pop, n = 3)

.   ID   CL   VC     KA KOUT IC50 FOO
. 1  1 1.05 47.8 0.8390 2.45 1.28   4
. 2  2 0.73 30.1 0.0684 2.51 1.84   6
. 3  3 2.82 23.8 0.1180 3.88 2.48   5

param(mod)

. 
.  Model parameters (N=3):
.  name value . name value
.  CL   1     | V    20   
.  KA   1     | .    .

What else can we do with idata?

Batches of simulations or sensitivity analyses

idata <- expand.idata(CL = seq(0.5, 1.5, 0.25))
idata

.   ID   CL
. 1  1 0.50
. 2  2 0.75
. 3  3 1.00
. 4  4 1.25
. 5  5 1.50

Note: this is the event + idata_set configuration

mod %>%
  idata_set(idata) %>%
  ev(amt = 100, ii = 24, addl = 2) %>%
  mrgsim(end = 120) %>% 
  plot(CP ~ ., logy=TRUE)

Your turn

• Description: sensitivity analysis in a PBPK model

• File name: exercises/sensitivity.R

data_set is the dosing equivalent to idata_set

data <- expand.ev(amt = c(100, 300, 1000), ii = 24, addl = 3)

head(data)

.   ID  amt ii addl evid cmt time
. 1  1  100 24    3    1   1    0
. 2  2  300 24    3    1   1    0
. 3  3 1000 24    3    1   1    0

Note: this is data_set configuration

mod %>%
  data_set(data) %>%
  mrgsim(end = 120) %>% plot(log(CP) ~ .)

data_set can also carry parameters

data <- expand.ev(
  amt = c(100,300,1000),
  ii = 24, addl = 3,
  CL = seq(0.5,1.5, 0.5)
)

head(data)

.   ID  amt ii addl  CL evid cmt time
. 1  1  100 24    3 0.5    1   1    0
. 2  2  300 24    3 0.5    1   1    0
. 3  3 1000 24    3 0.5    1   1    0
. 4  4  100 24    3 1.0    1   1    0
. 5  5  300 24    3 1.0    1   1    0
. 6  6 1000 24    3 1.0    1   1    0

carry_out

Copy from the input data to the output data

mod %>% carry_out(dose) %>% data_set(data, dose==300) %>% mrgsim()

. Model:  pk1 
. Dim:    2886 x 4 
. Time:   0 to 192 
. ID:     3 
.     ID time dose     CP
. 1:   2  0.0  300  0.000
. 2:   2  0.0  300  0.000
. 3:   2  0.2  300  2.712
. 4:   2  0.4  300  4.919
. 5:   2  0.6  300  6.712
. 6:   2  0.8  300  8.167
. 7:   2  1.0  300  9.345
. 8:   2  1.2  300 10.296

data_set and ev

head(data, n = 3)

.   ID  amt ii addl  CL evid cmt time dose
. 1  1  100 24    3 0.5    1   1    0  100
. 2  2  300 24    3 0.5    1   1    0  300
. 3  3 1000 24    3 0.5    1   1    0 1000

ev(amt = 100, ii = 24, addl = 3)

. Events:
.   time cmt amt evid ii addl
. 1    0   1 100    1 24    3

We have now seen 3 simulation setups

• mod + ev = ?
• mod + idata_set + ev = ?
• mod + data_set = ?

Wait a minute …

head(data)

.   ID  amt ii addl  CL evid cmt time dose
. 1  1  100 24    3 0.5    1   1    0  100
. 2  2  300 24    3 0.5    1   1    0  300
. 3  3 1000 24    3 0.5    1   1    0 1000
. 4  4  100 24    3 1.0    1   1    0  100
. 5  5  300 24    3 1.0    1   1    0  300
. 6  6 1000 24    3 1.0    1   1    0 1000

Generating data sets

Recall that data sets are just plain old data.frames in R. Feel free to code these up in any way that is convenient for you.

In our experience, the following helper functions cover many (not every) common needs for building the data sets.

• expand.ev
• as_data_set
• ev_rep
• ev_assign

expand.ev

• Like expand.grid … gives “all combinations”

expand.ev(ID = 1:2, amt = c(100,200))

.   ID amt evid cmt time
. 1  1 100    1   1    0
. 2  2 100    1   1    0
. 3  3 200    1   1    0
. 4  4 200    1   1    0

as_data_set

• Combine event objects into a single data set

data <- as_data_set(
  ev(amt = 100, ii = 12, addl = 19, ID = 1:2),
  ev(amt = 200, ii = 24, addl = 9,  ID = 1:3),
  ev(amt = 150, ii = 24, addl = 9,  ID = 1:4)
)

.   ID time cmt evid amt ii addl
. 1  1    0   1    1 100 12   19
. 2  2    0   1    1 100 12   19
. 3  3    0   1    1 200 24    9
. 4  4    0   1    1 200 24    9
. 5  5    0   1    1 200 24    9
. 6  6    0   1    1 150 24    9
. 7  7    0   1    1 150 24    9
. 8  8    0   1    1 150 24    9
. 9  9    0   1    1 150 24    9

ev_rep

• Replicate an event object for several IDs

e <- ev(amt = 100, ii = 12, addl = 14)
ev_rep(e, ID = seq(5))

.   time cmt amt evid ii addl ID
. 1    0   1 100    1 12   14  1
. 2    0   1 100    1 12   14  2
. 3    0   1 100    1 12   14  3
. 4    0   1 100    1 12   14  4
. 5    0   1 100    1 12   14  5

e <- ev(amt = 100, ii = 12, addl = 14, ID = seq(5))

ev_rep

e <- seq(ev(amt = 100), wait = 36, ev(amt = 50, ii = 24, addl = 4))
e

. Events:
.   time cmt amt evid ii addl
. 1    0   1 100    1  0    0
. 2   36   1  50    1 24    4

ev_rep(e, ID = seq(2))

.   time cmt amt evid ii addl ID
. 1    0   1 100    1  0    0  1
. 2   36   1  50    1 24    4  1
. 3    0   1 100    1  0    0  2
. 4   36   1  50    1 24    4  2

ev_assign

pop <- expand.idata(GROUP = c(1,2), ID = 1:3)
head(pop)

.   ID GROUP
. 1  1     1
. 2  2     2
. 3  3     1
. 4  4     2
. 5  5     1
. 6  6     2

e1 <- ev(amt = 100, ii = 24, addl = 9)
e2 <- mutate(e1, amt = 200)

data <- ev_assign(
  l = list(e1,e2), 
  idata = pop, 
  evgroup = "GROUP",
  join = TRUE
)

head(data)

.   time cmt amt evid ii addl ID GROUP
. 1    0   1 100    1 24    9  1     1
. 2    0   1 200    1 24    9  2     2
. 3    0   1 100    1 24    9  3     1
. 4    0   1 200    1 24    9  4     2
. 5    0   1 100    1 24    9  5     1
. 6    0   1 200    1 24    9  6     2

Your turn

• Description: simulate GCSF PK/PD profiles for weight-based dosing

• File name: exercises/gcsf.R

Create complex events - 1

What’s going to happen?

e1 <- ev(amt = 200) 

e2 <- ev(amt = 100, time = 24, ii = 24, addl = 2)

c(e1, e2)

Create complex events - 1

What’s going to happen?

Create complex events - 1

Combine two events

e1 <- ev(amt = 200) 

e2 <- ev(amt = 100, time = 24, ii = 24, addl = 2)

c(e1, e2)

. Events:
.   time cmt amt evid ii addl
. 1    0   1 200    1  0    0
. 2   24   1 100    1 24    2

Create complex events - 2

What’s going to happen?

e1 <- ev(amt = 200, ii = 12, addl = 2) 

e2 <- ev(amt = 100, ii = 24, addl = 4)

seq(e1, e2)

Create complex events - 2

What’s going to happen?

Create complex events - 2

Put two events in a sequence

e1 <- ev(amt = 200, ii = 12, addl = 2) 
e2 <- ev(amt = 100, ii = 24, addl = 4)

seq(e1, e2)

. Events:
.   time cmt amt evid ii addl
. 1    0   1 200    1 12    2
. 2   36   1 100    1 24    4

Create complex events - 3

What is going to happen?

e1 <- ev(amt = 200) 

e2 <- ev(amt = 100, ii = 24, addl = 2)

seq(e1, wait = 36, e2)

Create complex events - 3

What is going to happen?

Create complex events - 3

Wait before starting the next part of the regimen

e1 <- ev(amt = 200) 

e2 <- ev(amt = 100, ii = 24, addl = 2)

seq(e1, wait = 36, e2)

. Events:
.   time cmt amt evid ii addl
. 1    0   1 200    1  0    0
. 2   36   1 100    1 24    2

Your turn

• File name: exercises/z-pak.R

• Section name: Z-Pak

We’re stil working on this setup

model %>% intervention %>% simulate %>% take-a-look

model:

• Load a model with mread or mread_cache
• Use the internal library with modlib("<model-name>")
• Check model parameters with param(mod)
• Check model initial conditions with init(mod)
• View model code with see(mod)

intervention:

• ev(...): amt, cmt, time, ii, addl, rate
• Different ways to combine event objects

Your turn

• Description: Population PK of azithromycin

• File name: exercises/population_azithromycin.R

Your turn

• Description: simulate PK profiles using empirical Bayes estimates from a meropenem model fit in NONMEM

• File name: exercises/sim_ebe.R

Let’s make this model

• What are the parameters?
• What are the compartments?
• What are the differential equations?

Code up the model

Block $PARAM [R] declare and initialize model parameters

$PARAM CL = 1, VC = 20, KA=1.2

• Model “parameters” with name = value format
• Separate by , or <newline>
• Use any name except for words in mrgsolve:::reserved()
• Values are evaluated by R parser
• Multiple blocks allowed
• Use name anywhere in the model
• Updatable in R but read-only in the model specification file
  • Consider $FIXED if appropriate
• Most often, you will want to match names in $PARAM with the names in your input data sets

Block $CMT [txt] declare model compartments

• Specify name and number of compartments in the model
• Use any name except for those listed in mrgsolve:::reserved()
• See also or $INIT

Example $CMT

$CMT GUT CENT

Example $INIT

$INIT GUT = 0, CENT = 10

Block $ODE [C++] write differential equations

• Every compartment needs an equation
• Form:dxdt_CMT = ratein - rateout;
• Use CMT for compartment amounts, parameters, or other variables
• This block of code gets called repeatedly so be wise

$ODE
dxdt_GUT  = -KA*GUT;
dxdt_CENT =  KA*GUT - (CL/VC)*CENT;

Derive new variables in your model

• Your model code is written in C++, which is a language
• When you want a new variable, you must say what it is
• $MAIN, $ODE, $TABLE, $GLOBAL

To get a numeric variable, use double

double x = 5.4;

Other types you might use

bool cured = false;
int i = 2;

If in doubt, use double; it’s what you want most of the time.

Block $MAIN [C++] covariate model & initials

For example

$MAIN
double CL = TVCL*pow(WT/70,0.75);
double VC = TVVC*pow(0.86,SEX);
RESP_0 = KIN/KOUT;

In this example

• TVCL, KIN, KOUT, WT, SEX were declared in $PARAM
• There is a response compartment RESP

C++ expressions and functions

if(a == 2) b = 2;
if(b <= 2) {
  c=3;
} else {
  c=4;
}
double d = pow(base,exponent);
double e = exp(3);
double f = fabs(-4);
double g = sqrt(5);
double h = log(6);
double i = log10(7);
double j = floor(4.2);
double k = ceil(4.2);

Lots of help on the web http://en.cppreference.com/w/cpp/numeric/math/tgamma

Preprocessor directives (#define)

$GLOBAL
#define CP (CENT/VC)
#define g 9.8
$ODE
double INH = CP/(EC50+CP);
dxdt_RESP = KIN*(1-INH) - KOUT*RESP;
$CAPTURE CP

Introduction to POPULATION simulation

• Like NONMEM, use $OMEGA and $SIGMA to define covariance matrices for IIV and RUV, respectively

• Like NONMEM, use ETA(n) and EPS(m) for realized random effects drawn from $OMEGA and $SIGMA respectively
  • BUT …. use labeled ETAs instead

• Like NONMEM, mrgsolve recognizes ID as a subject indicator for simulating a new ETA

• Like NONMEM, mrgsolve allows you to write an error model as a function of EPS(m) and any other calculated value in your model

Block $OMEGA and $SIGMA [txt]

• Specify random effect variance/covariance matrices
  • \(\eta_i \sim N\left(0, \Omega \right)\)
  • \(\varepsilon_j \sim N \left(0, \Sigma \right)\)
• NONMEM-style input - lower triangle
• Multiple $OMEGA and $SIGMA are allowed and each may be named
• How you enter these blocks defines a that mrgsolve will enforce later on
• Diagonal matrix (3x3 in this example)

$OMEGA 0.1 0.2 0.3

More $OMEGA and $SIGMA

• Options are @ delimited & must be on different line from the matrix data
• Block matrix (2x2 in this example)

$OMEGA @block
0.1 0.0947 0.2

• Block matrix with correlations (assumes @block)

$OMEGA @cor
0.1 0.67 0.2

Block $OMEGA and $SIGMA

• Multiple $OMEGA and $SIGMA are allowed; each may be named

$OMEGA @name PK
0 0 0
$OMEGA @name PD
0 0

Users are encouraged to add labels

$OMEGA @name PK @labels ECL EVC EKA
0 0 0

About @ macros

• Used to indicate options associated with different code blocks
• Most options are for $OMEGA and $SIGMA
• @ should appear at the start of the line and the entire line is reserved for options
• A line may contain multiple options: @cor @name PK

Two forms:

• Boolean: @block means block=TRUE
• Name/value: @labels ECL EVC means labels=c("ECL", "EVC")

Closed-form PK models with $PKMODEL

• Set ncmt to 1 or 2
• Set depot to TRUE for extravascular dosing compartment
• Set trans to change the names of required parameters
• Use $CMT to declare 1 to 3 compartments as appropriate

$PKMODEL cmt = "GUT CENT PERIPH", depot=TRUE
$PARAM CL = 1 , V2 = 30, Q = 8, V3 = 400, KA=1.2

More $PKMODEL

$PKMODEL ncmt=1, depot=TRUE
$CMT GUT CENT
$PARAM TVCL = 1 , TVV = 30, KA=1.2
$OMEGA @labels ECL EVC
1 1
$MAIN
double CL = TVCL*exp(ECL);
double V  = TVV *exp(EVC);

Bioavailability, dosing lag time, and infusions

• To change the bioavability of doses administered to a compartment, set F_CMT in $MAIN

• To add a lagtime for doses administered to a compartment, set ALAG_CMT in $MAIN

• To set the duration of infusion, set D_CMT in $MAIN
  • Also, use rate = -2 in your data set or event object
• To set the infusion rate, set R_CMT in $MAIN
  • ALso, use rate = -1 in your data set or event object

Read in a model object with caching

First time to read

mod <- mread_cache("simple", "../../model")

Next time to read

mod <- mread_cache("simple", "../../model")

Model build directory (soloc)

It can be helpful to set

options(mrgsolve.soloc = "my_build_directory")

mod <- mread("simple", "model", soloc = "my_build_directory")

Different ways to call mread

mod <- mread("csa.txt", "../../model")

and

mod <- mread("../../model/csa.txt")

Inline model specification

We haven’t covered the specifics of coding a model yet

code <- '
$PARAM CL = 1, V = 20
$PKMODEL cmt = "CENT"
'

mod <- mcode("dont_do_this", code)

mod <- mcode_cache("seriously_dont", code)

Controlling output - request

mod %>% mrgsim() %>% head(n = 2)

.   ID time GUT CENT CP
. 1  1    0   0    0  0
. 2  1    1   0    0  0

mod %>% Req(CP) %>% mrgsim() %>% head(n = 2)

.   ID time CP
. 1  1    0  0
. 2  1    1  0

Controlling output - obsonly

mod %>% ev(amt=1) %>% mrgsim() %>% head(n = 2)

.   ID time GUT CENT CP
. 1  1    0   0    0  0
. 2  1    0   1    0  0

mod %>% ev(amt=1) %>% obsonly() %>% mrgsim() %>% head(n = 2)

.   ID time       GUT      CENT         CP
. 1  1    0 0.0000000 0.0000000 0.00000000
. 2  1    1 0.3011942 0.6782976 0.03391488

Sensitivity analysis

• Ad-hoc
  • Sequential: systematically vary the parameter value
  • Random: sample from a parametric or non-parametric distribution
• Local
  • small changes around a nominal value; one at a time
• Global
  • explore the whole parameter space; includes interactions

Local sensitivity analysis - sensFun()

We use the function sensFun from the package FME (FME::sensFun)

library(FME)

?sensFun

• To do local sensitivity analysis, we need
  • A function that accepts parameters as an argument and returns simulated data
  • Some parameters to test
  • An output variable of interest

sens <- sensFun(my_function, my_parameters, sensvar = my_variable_name)

• sensFun will take the parameter values we passed and manipulate them by small amounts to determine the local sensitivity

First, look at the result

A simulation scenario of interest

mod <- modlib("irm1", end = 72, delta = 0.1)

mod %>%
  ev(amt = 100) %>% 
  mrgsim(end = 72, delta = 0.1) %>% 
  plot

Create myfun to capture the scenario

• This is the same workflow, but we wrap it in a function and allow parameters (and the model) to be passed in as arguments

myfun <- function(pars, mod) {
  mod %>% 
    param(pars) %>% 
    ev(amt = 100) %>% 
    mrgsim(end = 72) %>%
    as.data.frame %>% 
    select(-ID)
}

• Also, we don’t plot the result, but rather return a data frame

Test out myfun

myfun(param(mod), mod) %>% head()

.   time       EV1      CENT     PERIPH     RESP EV2        CP
. 1  0.0   0.00000  0.000000 0.00000000 5.000000   0 0.0000000
. 2  0.0 100.00000  0.000000 0.00000000 5.000000   0 0.0000000
. 3  0.1  90.48374  9.444382 0.04780971 4.902817   0 0.4722191
. 4  0.2  81.87308 17.851273 0.18293635 4.688450   0 0.8925637
. 5  0.3  74.08182 25.323260 0.39389717 4.426064   0 1.2661630
. 6  0.4  67.03200 31.953048 0.67040259 4.150992   0 1.5976524

Pick some parameters

pars <- param(mod)[c("CL", "V2", "V3", "KA")]

pars

. $CL
. [1] 1
. 
. $V2
. [1] 20
. 
. $V3
. [1] 10
. 
. $KA
. [1] 1

Run the analysis

sens <- sensFun(myfun, pars, sensvar = "CP", mod = mod)

head(sens, n=3)

.     x var           CL         V2            V3        KA
. 1 0.0  CP  0.000000000  0.0000000  0.000000e+00 0.0000000
. 2 0.0  CP  0.000000000  0.0000000  0.000000e+00 0.0000000
. 3 0.1  CP -0.002544135 -0.9924278 -3.331475e-05 0.9506995

The default plot method

plot(sens, legpos = "right")

Re-shape sensitivity analysis results

sens %>% 
  gather(variable,value,CL:KA) %>% 
  ggplot(aes(x,value,col=variable)) + geom_line(lwd=1) + theme_bw()

Global senstivity analysis

• All the parameters vary at once
• Investigate the whole parameter space
  • Or reasonable parameter space
• Evaluate both individual parameters as well as interactions between parameters
• Much greater computational burden

We’ll look at Sobol’s method

Global senstivity analysis

• Sobol Sensitivity Analysis: A Tool to Guide the Development and Evaluation of Systems Pharmacology Models.
  • Zhang XY, Trame MN, Lesko LJ, Schmidt S.
  • CPT Pharmacometrics Syst Pharmacol. 2015 Feb;4(2):69-79. doi: 10.1002/psp4.6.
  • PMID: 27548289
  • https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5006244/

Global senstivity analysis

• Sobol Sensitivity Analysis: A Tool to Guide the Development and Evaluation of Systems Pharmacology PBPK Models.
  • Zhang XY, Trame MN, Lesko LJ, Schmidt S.
  • CPT Pharmacometrics Syst Pharmacol. 2015 Feb;4(2):69-79. doi: 10.1002/psp4.6.
  • PMID: 27548289
  • https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5006244/

Sobol

• Variabillity in the output is decomposed, relating to variability in different inputs or combinations of inputs

Case study: ciprofloxacin PBPK

• fluoroquinolone antibiotic
  • UTI , uncomplicated cystitis, SSTI, respiratory tract infection
  • Pseudomonas aeroginosa
  • Bacillus anthracis
• PBPK model
  • A whole-body physiologically based pharmacokinetic (WB-PBPK) model of ciprofloxacin
  • https://www.ncbi.nlm.nih.gov/pubmed/27578330
  • Code in ./model/cipro.cpp
• AUC:MIC ratio
  • The ciprofloxacin target AUC : MIC ratio is not reached in hospitalized patients with the recommended dosing regimens
  • https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3555057

Cipro PK: warm-up

mod <- mread("cipro", "../../model")

• What are the parameters?
• What are the compartments?
• What is the model code?

Cipro PK: infusion doses

• Simulate 400 mg as IV infusion over 1 hour

• Plot cipro concentration over 12 hours

Cipro PK: infusion doses

• Simulate 400 mg as IV infusion over 1 hour

• Plot cipro concentration over 12 hours

• Simulate 400 mg one day of a q12h dosing

• Plot over 24 hours of dosing

Cipro PK: covariates

Explore influence of patient weight, sex, and creatinine clearance on cipro PK

• Weight: 40 to 140 kg
• Creatinine clearance: 90, 60, 45, 30, 15 mL/min
• Sex: male / female
• Inputs: idata_set, event object
• Outputs: ciprofloxacin concentration, covariate values

Cipro PK: AUC

From the Haeseker paper:

Mean (+/- SD) age was 66 (+/-17) years, the mean clearance corrected for bodyweight was 0.24 l h-1 kg^1 and the mean AUC was 49 mg l^1 h.

Table 2 shows CrCL, CL, V, and AUC by total daily ciprofloxacin dose

Let’s confirm the AUC result for each of the covariate groups in the previous step.

Use PKPDmisc::auc_partial to calculate AUC.

Cipro PK: organ-specific exposure

Profiles of unbound ciprofloxacin concentration in different tissues were predicted from the developed WB-PBPK model (Fig. 3). Kidney and lung were predicted to have higher exposures as compared to the other organs while muscle, brain and adipose were predicted to have relatively low exposures to ciprofloxacin.

Let’s look at exposure in different organs

Use ./model/cipro_conc.cpp model

Cipro PK: simulation with uncertainty

ciprofloxacin 400 mg IV over 2 hours

What is the distribution of AUC after the last dose on day 3 of q12h dosing?

Use ./model/cipro.cpp

Use ./data/cipro_post.RDS

Cipro PK: AUC:MIC

Simulate the distribution of 24-hour AUC:MIC for

• MIC: 0.5
  • 200 mg IV q12h over 1 hour with CRCL = 30 ml/min
  • 400 mg IV q12h over 1 hour, normal renal function
• MIC: 1
  • 400 mg IV q12h over 1 hour
  • 400 mg IV q8h over 1 hour

Check progress on the second day of dosing